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Abstract 

We  present  a  number  of  complete  sets  of  series  expansion  formulas  for  neutrino  oscillation 
probabilities  in  matter  of  constant  density  for  three  flavors.  In  particular,  we  study  expansions 
in  the  mass  hierarchy  parameter  a  =  and  mixing  parameter  S13  =  sin  $13  up  to 

second  order  and  expansions  only  in  a  and  only  in  S13  up  to  first  order.  For  each  type  of 
expansion  we  also  present  the  corresponding  formulas  for  neutrino  oscillations  in  vacuum.  We 
perform  a  detailed  analysis  of  the  accuracy  of  the  different  sets  of  series  expansion  formulas 
and  investigate  which  type  of  expansion  is  most  accurate  in  different  regions  of  the  parameter 
space  spanned  by  the  neutrino  energy  E,  the  baseline  length  L ,  and  the  expansion  parameters 
a  and  S13.  We  also  present  the  formulas  for  series  expansions  in  a  and  in  S13  up  to  first  order 
for  the  case  of  arbitrary  matter  density  profiles.  Furthermore,  it  is  shown  that  in  general  all 
the  18  neutrino  and  antineutrino  oscillation  probabilities  can  be  expressed  through  just  two 
independent  probabilities. 
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1  Introduction 


The  discovery  of  neutrino  oscillations  in  atmospheric,  solar,  and  reactor  neutrino  experiments 
has  turned  neutrino  physics  into  one  of  the  most  exciting  and  active  fields  of  particle  physics. 
By  now,  a  significant  amount  of  information  on  neutrino  properties  has  been  obtained.  The 
results  of  the  atmospheric  neutrino  experiments  [1-3]  and  the  K2K  accelerator  neutrino  ex¬ 
periment  [4]  have  allowed  the  determination  of  the  fundamental  neutrino  parameters  A 
and  d23  to  an  accuracy  of  about  30%  and  15%,  respectively,  while  the  solar  neutrino  experi¬ 
ments  [5]  and  the  KamLAND  reactor  neutrino  experiment  [6]  have  measured  the  parameters 
Am],  and  9 u  to  an  accuracy  of  about  15%.  In  fact,  this  means  that  neutrino  physics  is 
now  entering  an  era  of  precision  measurements  of  the  neutrino  oscillation  parameters.  Future 
experiments  with  superbeams  and  neutrino  factories  will  determine  the  “atmospheric”  and 
“solar”  neutrino  oscillation  parameters  to  an  accuracy  of  the  order  of  1%.  They  are  also 
expected  to  measure  the  elusive  leptonic  mixing  angle  6 13  for  which  at  present  only  an  upper 
limit  exists  [7],  or  to  put  a  more  stringent  limit  on  this  angle,  as  well  as  to  clarify  the  issues 
of  the  neutrino  mass  hierarchy  and  possibly  of  leptonic  CP  violation.  Important  information 
on  neutrino  oscillation  parameters  can  also  be  obtained  from  the  future  atmospheric,  solar, 
reactor,  and  supernova  neutrino  experiments. 

Increasing  accuracy  and  reach  of  the  present  and  especially  of  forthcoming  experiments 
put  forward  new  and  more  challenging  demands  to  the  theoretical  description  of  neutrino 
oscillations.  In  order  to  be  able  to  determine  the  fundamental  neutrino  parameters  from  the 
data  with  high  precision,  one  needs,  among  other  things,  very  accurate  theoretical  expressions 
for  the  probabilities  of  neutrino  oscillations  in  matter  and  in  vacuum.  While  these  proba¬ 
bilities  can,  in  principle,  be  calculated  numerically  with  any  requisite  accuracy,  it  is  highly 
desirable  to  have  also  analytic  expressions  for  them.  Such  analytic  expressions  would  reveal 
the  basic  dependence  of  the  neutrino  oscillation  probabilities  on  the  fundamental  neutrino 
parameters  and  on  the  characteristics  of  the  experiment  and  thus  facilitate  the  choice  of  the 
experimental  setup  as  well  as  the  analysis  of  the  data.  They  would  also  help  to  understand 
the  physics  underlying  various  flavor  transitions  and  to  resolve  the  parameter  degeneracies 
and  other  ambiguities,  such  as  fundamental  versus  matter-induced  CP  violation. 

The  purpose  of  this  paper  is  to  present  a  collection  of  approximate  analytic  formulas 
for  the  neutrino  oscillation  probabilities.  We  derive  a  number  of  complete  sets  of  series 
expansion  formulas  for  the  three-flavor  neutrino  oscillation  probabilities  in  matter  and  in 
vacuum.  The  probabilities  are  expanded  in  the  mass  hierarchy  parameter  a  = 
the  mixing  parameter  s13  =  sin  6*13,  or  in  both  of  them.  We  also  study  in  detail  the  accuracy 
of  the  obtained  expressions  in  different  regions  of  the  parameter  space  and  identify  the  “best 
choice”  in  each  case  of  interest. 

Before  proceeding  to  present  our  results,  we  give  here  a  brief  overview  of  the  previous 
work  on  the  subject.  Analytic  formulas  for  three-flavor  neutrino  oscillation  probabilities  have 
been  derived  in  a  number  of  papers.  Exact  formulas  for  the  neutrino  oscillation  probabilities 
in  vacuum  can  be  found,  e.g.,  in  Ref.  [8].  Exact  expressions  can  also  be  obtained  in  the 
case  of  three-flavor  neutrino  oscillations  in  matter  of  constant  density  [9-16].  However,  the 
corresponding  formulas  are  rather  complicated  and  not  easily  tractable.  This  also  applies 
to  the  exact  analytic  three-flavor  formulas  obtained  for  some  special  cases  of  non-uniform 
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matter  density:  linear  matter  density  [17]  and  exponentially  varying  matter  density  [18]. 

Approximate  analytic  formulas  for  three-flavor  neutrino  oscillation  probabilities  in  mat¬ 
ter  have  been  derived  in  a  number  of  papers,  see  e.g .,  Refs.  [19-24],  In  Refs.  [19,20],  the 
Mikheyev-Smirnov-Wolfenstein  (MSW)  resonances  [25]  in  matter  of  varying  density  were 
studied  assuming  large  separation  and  independence  of  the  high-density  and  low-density  res¬ 
onances.  In  Refs.  [20,21,23],  the  adiabatic  approximation  was  used  to  derive  three-flavor  neu¬ 
trino  oscillation  probabilities  in  matter  of  varying  density,  whereas  in  Ref.  [24],  the  Magnus 
expansion  for  the  time  evolution  operator  was  used,  which  is  equivalent  to  the  anti-adiabatic 
approximation. 

In  a  number  of  papers  an  approach  similar  to  ours  was  adopted,  i.e.,  the  neutrino  os¬ 
cillation  probabilities  were  expanded  in  the  small  parameters  a,  .S13,  or  in  both  of  them. 
In  Ref.  [26],  exact  analytic  expressions  for  three-flavor  neutrino  oscillation  probabilities  in 
matter  with  an  arbitrary  density  profile  were  obtained  in  the  limit  a  — >  0  by  reducing  the 
problem  to  an  effective  two-flavor  one.  In  Ref.  [27],  a  similar  approach  was  employed  to  ob¬ 
tain  the  neutrino  oscillation  probabilities  in  matter  of  arbitrarily  varying  density  in  the  limit 
$13  — >  0,  whereas  in  Refs.  [28,29],  expressions  up  to  first  order  in  S13  were  derived.  For  the 
case  of  matter  of  constant  density,  the  limit  a  — >  0  was  considered  in  Ref.  [30].  Expansions 
up  to  first  order  in  a  were  carried  out  in  Refs.  [31,32],  In  Refs.  [33,34],  both  the  solar  mass 
squared  difference  and  matter  effects  were  treated  as  perturbations  and  the  transition 

probabilities  up  to  first  order  in  them  were  derived.  Expansions  in  both  a  and  S13  were  used 
in  Refs.  [35-38].  We  note  that  the  neutrino  oscillation  probabilities  derived  in  most  of  the 
above-mentioned  papers  either  did  not  constitute  a  complete  set  (i.e.,  a  set  of  probabilities 
from  which  the  probabilities  in  all  channels  can  be  obtained),  or  contained  expressions  from 
which  some  terms  were  missing,  especially  in  the  case  of  the  probabilities  P(V/t  — >  ///t)  and 
P(z/M  — >  uT).  Thus,  to  the  best  of  our  knowledge,  our  study  is  the  first  one  in  which,  for  the 
case  of  matter  of  constant  density,  complete  and  consistent  expansions  in  a  and  Si3  to  second 
order  and  expansions  only  in  a  and  only  in  S13  up  to  first  order  are  performed. 

The  paper  is  organized  as  follows.  In  Sec.  2  we  set  the  general  formalism  and  notation, 
discuss  some  general  relations  satisfied  by  the  oscillation  probabilities,  and  show  that  all 
the  18  neutrino  and  antineutrino  oscillation  probabilities  can  be  expressed  through  just  two 
independent  probabilities.  In  Sec.  3  we  present  a  series  expansion  of  the  neutrino  oscillation 
probabilities  up  to  second  order  in  both  a  and  S13  (the  so-called  “double  expansion”),  whereas 
in  Secs.  4  and  5  we  consider  the  probabilities  expanded  up  to  first  order  in  a  and  S13, 
respectively  (“single  expansions”).  In  all  three  cases,  we  give  the  probabilities  for  matter  of 
constant  density  (Secs.  3.1,  4.1,  and  5.1)  and  in  vacuum  (Secs.  3.2,  4.2,  and  5.2).  We  also 
compare  our  formulas  with  the  corresponding  expressions  existing  in  the  literature  (when 
available),  pointing  out  agreements  and  disagreements.  In  Sec.  6,  we  discuss  the  qualitative 
behavior  of  the  neutrino  oscillation  probabilities,  the  relevance  of  matter  effects,  and  give  a 
detailed  evaluation  of  the  accuracy  of  the  various  formulas.  Furthermore,  we  comment  on  the 
application  of  our  probability  formulas  to  neutrino  oscillation  experiments.  We  summarize 
our  results  in  Sec.  7.  Finally,  several  methods  that  have  been  used  to  derive  the  formulas  are 
presented  in  the  appendices.  In  Appendix  A,  we  describe  the  perturbative  diagonalization 
of  the  effective  Hamiltonian  of  the  neutrino  system  in  the  case  of  matter  of  constant  density, 
while  in  Appendix  B,  the  details  of  the  perturbative  expansion  of  the  neutrino  evolution 


2 


equation  in  the  case  of  arbitrary  matter  density  profiles  are  given. 

2  General  formalism  and  notation 

We  consider  three-flavor  neutrino  oscillations  and  adopt  the  standard  parameterization  of  the 
leptonic  mixing  matrix  U  [39]: 

u  =  o23u5oV3ulo12 

(  C12C13  S12C13  <si3e_ll5cp\ 

=  -S\2C23  -  Ci2Sl3S23el5cp  C12C23  -  Sl2-Sl3S23el5cp  Ci3S23  •  (1) 

V  S12S23  -  Ci2Sl3C23e1,5cp  -C12S23  -  Sl2<Sl3C23e1,5cp  C13C23  ) 

Here  OtJ  is  the  orthogonal  rotation  matrix  in  the  ij-plane  which  depends  on  the  mixing 
angle  0lv  Us  =  diag(l,  1,  el5cp),  $cp  being  the  Dirac-type  CP- violating  phase,  =  sin 
and  Cij  =  cos  6 ij.  In  the  three-flavor  case  there  are  also,  in  general,  two  Majorana-type  CP- 
violating  phases;  however,  these  phases  do  not  affect  neutrino  oscillations,  and  therefore  will 
not  be  considered  here.  Without  loss  of  generality,  one  can  assume  all  the  mixing  angles  to 
lie  in  the  first  quadrant  (he.,  between  0  and  7t/2),  while  the  CP-violating  phase  5cp  is  allowed 
to  lie  in  the  interval  [0,  27t]  . 

Let  us  denote  by  Pap  =  P(ua  — >  up)  the  probability  of  transition  from  a  neutrino  flavor 
a  to  a  neutrino  flavor  /?,  and  similarly  for  antineutrino  flavors,  he.,  P^p  =  P{pQ  — >  Up). 
In  general,  the  three-flavor  neutrino  oscillation  probabilities  in  matter  Pap  depend  on  eight 
parameters  and  one  function: 

Pap  =  PapiAmli,  Aml1,912,913,923,SCp-,E,L,V(x)),  a,(3  =  e,fjL,T .  (2) 

Here  A m|-  =  m]  —  rrij  are  the  neutrino  mass  squared  differences,  E  is  the  neutrino  en¬ 
ergy,  L  is  the  baseline  length,  and  V (. x )  is  the  matter- induced  effective  potential,  x  e  [0,  L\ 
being  the  coordinate  along  the  neutrino  path.  The  neutrino  mass  squared  differences,  the 
leptonic  mixing  angles,  and  the  CP-violating  phase  are  fundamental  parameters  and  thus 
experiment- independent,  whereas  the  neutrino  energy,  the  baseline  length,  and  the  matter- 
induced  effective  potential  vary  from  experiment  to  experiment.  The  present  best-fit  values 
and  3<v  allowed  ranges  of  the  fundamental  neutrino  parameters  found  in  a  recent  global  fit 
of  the  neutrino  oscillation  data  [40]  are  summarized  in  Table  1.  Unless  otherwise  stated, 
all  calculations  in  the  present  paper  are  performed  for  the  following  values  of  the  neutrino 
parameters:  Am|i  =  7  •  10”5eV2,  0\2  =  33°,  and  023  =  45°.  For  the  atmospheric  mass 
squared  difference  Am] ,  we  adopt  the  current  best-fit  value  given  by  the  Super-Kamiokande 
Collaboration,  IAtttJ-J  =  2  •  10”3  eV2  [2],  which  is  slightly  smaller  than  the  value  given  in  Ta¬ 
ble  1.  The  sign  of  Am]j  is  related  to  the  neutrino  mass  hierarchy:  for  the  normal  (inverted) 
hierarchy  one  has  Am]j  >  0  (Am]j  <  0).  For  the  leptonic  mixing  angle  633  we  allow  values 
below  the  90%  C.L.  upper  bound  found  in  the  global  fit  of  the  neutrino  oscillation  data  [40] 
for  \Amll\  fixed  at  2  •  10_3eV2: 

6*13  <  10.8°  ,  or  si3  =  sin  0,3  <  0.19  ,  or  sin2  20, 3  <  0.14  .  (3) 

For  the  CP  violation  phase  $cp  we  allow  values  between  0  and  27t. 
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Parameter 

Best-fit  value 

Range  (3<r) 

A  m221 

6.9-  10"5  eV2 

(5.4  4-  9.5)  •  10~5  eV2 

|An4| 

2.6  •  10-3  eV2 

(1.5  -T-  3.7)  •  10-3  eV2 

^12 

33.2° 

28.6°  -T-  38.6° 

^13 

4.4° 

0^  13.4° 

QA 

to 

CO 

46.1° 

33.8°  4- 58.1° 

<^CP 

- 

0  4-2tt 

Table  1:  Present  best-fit  values  and  3a  allowed  ranges  of  the  fundamental  neutrino  param¬ 
eters  from  a  three- flavor  fit  to  global  neutrino  oscillation  data  [40]. 


Inspecting  the  values  of  the  fundamental  neutrino  parameters  in  Table  1,  one  can  identify 
two  natural  candidates  for  small  expansion  parameters  of  the  neutrino  oscillation  probabil¬ 
ities.  These  are  the  small  lcptonic  mixing  angle  9Vi  (or,  equivalently,  S13)  and  the  mass 
hierarchy  parameter 


Am |i 

A  m§! 


~  0.026 , 


(0.018)  0.021  <  a  <  0.036  (0.053)  at  90%  C.L.  (3<r) , 


(4) 


where  we  have  taken  the  current  best-fit  value  and  allowed  ranges  from  Ref.  [40].  In  the 
following,  we  will  derive  a  number  of  formulas  for  series  expansions  of  the  neutrino  oscillation 
probabilities  in  these  small  quantities.  Comparing  Eqs.  (3)  and  (4),  one  realizes  that  current 
data  constrain  the  parameter  a  to  a  relatively  narrow  range,  while  .S13  is  only  bounded  from 
above  and  might  be  significantly  larger  or  smaller  than  a.  The  relative  size  of  these  two 
expansion  parameters  will  be  important  for  the  validity  of  a  given  type  of  expansion. 

In  order  to  find  the  neutrino  oscillation  probabilities  for  a  given  experimental  setup,  one 
has,  in  general,  to  solve  the  Schrodinger  equation  for  the  neutrino  vector  of  state  in  the  flavor 
basis  | u(t))  =  ( ue(t )  v^t)  z/T(f))r: 


with  the  effective  Hamiltonian 


u{t))  =  H\u{t )) 


(5) 


H  ~ 


1 

2 E 


U  diag(0,  A m221,  A m^JlA  +  diag(V,  0,  0) . 


(6) 


Here  V  is  the  charged-current  contribution  to  the  matter-induced  effective  potential  of  ue  [25]. 
We  have  disregarded  the  neutral-current  contributions  to  the  neutrino  potentials  in  matter, 
since  they  are  the  same  for  uei  and  vT  1  and  so  do  not  affect  neutrino  oscillations.  Note  that 
Eq.  (6)  holds  for  neutrinos,  whereas  for  antineutrinos  one  has  to  perform  the  replacements 


U  C*  ,  V  — >  —V  . 


(7) 


to  tiny  radiative  corrections  [41]  which  are  negligible  except  at  very  high  densities  available,  e.g., 
inside  supernovae. 
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The  potential  V  (x)  is  given  in  convenient  units  by 

V(x)  ~  7.56  x  10~14  ( YJx)  eV  ,  (8) 

\g/cm 6  ) 

where  p(x)  is  the  matter  density  along  the  neutrino  path  and  Ye(x)  is  the  number  of  electrons 
per  nucleon.  For  the  matter  of  the  Earth  one  has,  to  a  very  good  accuracy,  Ye  ~  0.5. 

For  many  practical  applications  (such  as  long-baseline  accelerator  experiments,  as  well 
as  oscillations  of  atmospheric,  solar,  and  supernova  neutrinos  inside  the  Earth  when  they 
do  not  cross  the  Earth’s  core)  it  is  a  very  good  approximation  to  assume  that  the  matter 
density  along  the  neutrino  trajectory  is  constant  (see,  e.g.,  Refs.  [42-44]).  Typical  values  for 
the  matter  density  are  pcmst  —  3g/cm3  in  the  Earth’s  crust  and  pmantie  —  4.5g/cm3  in  its 
mantle.  In  situations  where  the  neutrinos  also  cross  the  Earth’s  core  or  for  strongly  varying 
matter  density  profiles  like  those  inside  the  Sun  or  supernovae,  the  constant  matter  density 
approximation  is  not  valid. 

The  neutrino  oscillation  probabilities  can  be  found  as  Pa/3  =  \Spa(t,  to)  |2,  where  S(t,t0 ) 
is  the  evolution  matrix  such  that 

\u(t))  =  S(t,t0)\u(t0)} ,  S(t0,t0)  =  l.  (9) 

Note  that  S(t,t0 )  satisfies  the  same  Schrodinger  equation,  Eq.  (5),  as  | u(t)).  In  the  case 
of  matter  of  constant  density,  the  evolution  matrix  can  be  obtained  by  diagonalizing  the 
Hamiltonian  in  Eq.  (6)  according  to  H  =  U'HU'^ ,  where  U'  is  the  leptonic  mixing  matrix  in 
matter  and  H  =  diag(E1,  E2,  E:i) .  The  evolution  matrix  is  then  given  by 

3 

Spa(t,t0)  =  ^(U'aiyU,l3ie~lEiL ,  a,/3  =  e,/i,r,  (10) 

i=  1 

where  we  have  identified  L  =  t  —  t0. 

Before  presenting  our  results  for  the  neutrino  oscillation  probabilities,  we  discuss  some  of 
their  general  properties  as  well  as  relations  between  them.  First,  we  note  that  Eq.  (7)  implies 
that  one  can  relate  the  oscillation  probabilities  for  antineutrinos  to  those  for  neutrinos  by 

Pap  =  Pap(Si cp -£cp,  R -R) ,  a,  (3  =  e,  n,r .  (11) 

Second,  in  general,  i.e.,  both  in  vacuum  and  in  matter  with  an  arbitrary  density  profile,  it 
follows  from  the  unitarity  of  S(t,to)  (conservation  of  probability)  that 

^  '  Paj3  ^  ^  Pap  1  ,  CX,  f3  e,  /i,  T  .  (12) 

a  /3 

These  relations  imply  that  five  out  of  the  nine  neutrino  oscillation  probabilities  can  be  ex¬ 
pressed  in  terms  of  the  other  four  [45]. 

Besides  these  general  properties,  there  exists  an  additional  symmetry  due  to  the  specific 
parameterization  of  the  mixing  matrix  given  in  Eq.  (1)  and  the  fact  that  the  rotation  matrix 
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O23  commutes  with  the  matter  potential  term  of  the  Hamiltonian  in  Eq.  (6).  It  is  easy  to 
show  that  the  evolution  matrix  can  be  written  as 


S ( t ,  t0)  =  023S'(t,  t0)Oj3  (13) 

where  S'(t,t0 )  does  not  depend  on  6*23.  This  can  be  used  to  prove  some  useful  relations 
between  the  probabilities.  Let  us  denote  2 

Pa/3  =  Panels  C23  ,  sin  2023  -> -sin  2023) ,  a, /?  =  e, /1,  r  .  (14) 

Using  Eqs.  (13)  and  (14),  one  can  readily  show  that 

Per  =  Pen  ,  Prn  =  P^T  ,  Ptt  =  P nn  ,  (15) 


while  Pee  turns  out  to  be  independent  of  6*23. 

Out  of  the  three  conditions  in  Eq.  (15),  only  two  are  independent,  as  each  of  them 
can  be  derived  from  the  other  two  and  the  unitarity  conditions  (12).  Hence,  the  number 
of  independent  neutrino  oscillation  probabilities  is  reduced  to  two.  Thus,  we  come  to  the 
important  conclusion  that  all  the  nine  neutrino  oscillation  probabilities  can  be  expressed 
through  just  two  independent  probabilities  provided  that  their  dependence  on  the  mixing 
angle  623  is  known.  However,  the  choice  of  these  independent  probabilities  is  restricted: 
they  should  not  include  Pee,  which  is  independent  of  6*23;  nor  should  they  be  a  pair  of  the 
probabilities,  which  are  T-reverse  of  each  other,  or  go  into  each  other  (or  T-reverse  of  each 
other)  under  the  transformation  s|3  <-»■  c|3  ,  sin  2023  — 1 ►  —  sin  2023- 

One  possible  choice,  and  the  one  that  we  will  use,  is  Pe/1  and  P^T.  For  completeness,  we 
give  here  the  expressions  for  all  the  other  neutrino  oscillation  probabilities  in  terms  of  these 
two.  Using  Eqs.  (12)  and  (15),  one  easily  finds 


p 

1  ee 

—  1  —  (Pen  +  Pen)  ) 

(16) 

P 

1  er 

=  Pen  > 

(17) 

P 

1  ne 

u 

•of 

+ 

of 

1 

(IS) 

p 

1  w 

=  1  -  p  -  p 

-1-  J  efl  1  [AT  ) 

(19) 

p 

1  re 

=  P  -L  P  -  p 

J  e/Z  U  1  fAT  1  /AT  5 

(20) 

P 

1  Tfl 

=  Pnr  > 

(21) 

P 

1  TT 

=  1  —  ( Pen  +  Pnr)  ■ 

(22) 

In  addition  to  the  above  relations,  one  can  study  the  transformations  of  the  neutrino 
oscillation  probabilities  under  the  time  reversal  Pap  — >  Pga .  It  can  be  shown  [28]  that  in 
matter  with  an  arbitrary  density  profile 

Pf3a  =  Pag(SCp  ->  -6cp,  V(x)  — >  Vrev(x)) ,  a,  (3  =  e,  n,  r ,  (23) 

2The  transformation  in  Eq.  (14)  can  be  achieved,  e.g.,  through  the  substitution  023  — >  023  +  tt/2  or 
023  — ►  023  +  37r/2. 
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where  Viev(x)  is  the  “reverse”  potential,  which  corresponds  to  the  interchanged  positions 
of  the  neutrino  source  and  the  detector.  In  the  case  of  symmetric  matter  density  profiles 
(including  matter  of  constant  density),  VTev(x)  =  V(x),  and  Eq.  (23)  simplifies  to 

P/3a  =  Paf 3(^CP  “ >  — <5cp)  ,  «, /?  =  e, /i,  T  .  (24) 

While  Eq.  (24)  does  not  further  reduce  the  number  of  independent  probabilities,  it  yields 
relations  which  can  be  useful  for  cross-checking  the  formulas  for  Pap  in  the  case  of  matter 
with  symmetric  density  profiles. 

By  applying  the  rule  given  in  Eq.  (11),  one  can  obtain  from  Eqs.  (16)-(22)  the  cor¬ 
responding  probabilities  for  the  antineutrino  oscillations.  Thus,  the  expressions  for  all  18 
probabilities  of  neutrino  and  antineutrino  oscillations  can  be  found  from  the  formulas  for 
just  two  independent  neutrino  oscillation  probabilities,  which,  as  was  already  mentioned,  we 
choose  to  be  Pefl  and  P^r.  In  order  to  be  more  explicit,  we  will  in  some  cases  also  give 
formulas  for  additional  neutrino  oscillation  channels. 

In  the  following  Secs.  3,  4,  and  5,  we  give  our  results  for  the  various  series  expansions  of 
the  neutrino  oscillation  probabilities.  We  will  adopt  the  following  abbreviations: 

(25) 

(26) 

3  Series  expansion  up  to  second  order  in  a  and  Si3 


.  A moi  L 

A  =  31 


A  = 


AE  ’ 

2  EV  VL 
2A 


Amli 


3.1  Matter  of  constant  density 


In  this  section,  we  present  the  series  expansion  formulas  for  three-flavor  neutrino  oscilla¬ 
tion  probabilities  in  matter  of  constant  density  up  to  second  order  in  both  a  and  .s'13.  The 
probabilities  are  calculated  by  diagonalizing  the  Hamiltonian  (6)  up  to  second  order  in  these 
parameters,  as  described  in  Appendix  A.l.  We  find  for  the  eigenvalues  of  the  Hamiltonian 


Ei  ~ 

E2  ~ 
E3  ~ 


Arr4i 

2  E 

2  E 

2  E 


4  T  a  T  s 


A  ^ 2 


12^°13A-1 
2 


T  ot 


sin2  2012 
4  A 


CX  C^2  CX 


sin  2012 


4  A 


1  —  s 


A 


13 


A  -  1 


(27) 

(28) 
(29) 


Calculating  the  mixing  matrix  in  matter  U'  as  described  in  Appendix  A.l  and  using  Eq.  (10) 
for  the  evolution  matrix  S,  it  is  straightforward  (although  somewhat  tedious)  to  derive  the 
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following  expressions  for  the  neutrino  oscillation  probabilities: 


Pp e  —  1  —  a  sin  20 


2  2  sin2  ,  2  sin2(A  -  1) A 


-4s 


(A  —  l)2  ’ 


n  2  •  2  nn  2  Sill2  AA  2  2  sill2(A-l)A 

Pen  =  oc  Sin  2 012  C23  +  4  s13  s23  ^  _  ^2 


2  a  Si3  sin  2012  sin  2  023  cos  (A  —  $cp) 


sinAA  sin  (A  —  1)A 
A  A- 1 


sin2  AA  2  2  sin2(A  —  1)A 


PeT  =  a  sin  2012  s23 - — —  +  4  sj3  c23 


2  a  S13  sin  20i2  sin  2 023  cos(A  —  5cp) 


(A -1)2 

sin  AA  sin(A  —  1)A 
A  A^l 


Paa  =  1  —  sin2  2023  sin2  A  +  a  c?9  sin2  2023  A  sin  2  A 


sin2  AA 


—  a  sin  2012  c23 - — - a  cj2  sin  2 023  A  cos  2  A 


1  9  9  9  „  /  sin  A  A  ,  ,  A 

H - —  a2  sin2  20i2  sin“  2023  I  sin  A  — - —  cos(A  —  1)  A - sin  2A 

2A  \  A  2 


_  2  2  sin2  (A  —  1)A 

4513523  (A -l)2 

- —  s?3  sin2  2023  f  sin  A  cos  AA  S*n^ — —  —  A  sin  2 A 

A  -  1  13  3  \  A  - 1  2 


2  a  si3  sin  20i2  sin  2 023  cos  5cp  cos  A 


sin  AA  sin(A  —  1)A 
“ 4  A-  1 


2  f  sm  A  /\  \ 

+  — — -asi3  sin2012  sin2023  cos2023  cos^cp  sin  A  (A  sin  A - — —  cos(A  —  1)AJ  , 


PUT  =  sin2  2  023  sin2  A  —  a  c?9  sin2  2023  A  sin  2  A  +  a2  ct,  sin2  2  023  A2  cos  2  A 


- -  a2  sin2  2019  sin2  2 023  (  sin  A 


sin  A  A  i  A  A  . 

— - —  cos(A  —  1)A - sm2A 

A  2 


2  9  9  f  sin  (A  —  1)A  A  . 

+  — - sf.>  sm2  2023  (  sm  A  cos  AA - - - A  sm  2 A 

A  —  1  13  3  V  A  —  1  2 


2ccsi3  sin2012  sin2023  sin  sin  A 


sin  AA  sin(A  —  1)A 
A  A- 1 


— — ^-a<si3  sin  2012  sin  2023  cos2023  cos^cp  sin  A  (  A  sin  A 


sin  A  A 


(34) 


Formally,  our  calculations  are  based  upon  the  approximations  a,  s  t 3  -C  1  and  no  explicit 
assumptions  about  the  values  of  L/ E  are  made.  However,  we  remark  that  the  series  expansion 
formulas  (30)-(34)  are  no  longer  valid  as  soon  as  aA  =  Am^L  /  (4E)  becomes  of  order  unity, 
i.e.,  when  the  oscillatory  behavior  due  to  the  “solar”  mass  squared  difference  Arn^  becomes 
relevant.  This  can  happen  for  very  long  baselines  and/or  very  low  energies.  See  also  Sec.  6 
for  a  detailed  discussion  of  the  accuracy  of  these  formulas. 

From  Eqs.  (27)-(29)  one  can  see  that  in  vacuum  ( A  =  0)  and  at  the  atmospheric  reso¬ 
nance  ( A  =  1)  the  expressions  for  the  eigenvalues  are  divergent,  and  one  would  expect  the 
expansion  to  break  down.  In  these  cases,  two  out  of  the  three  eigenvalues  of  the  unperturbed 
Hamiltonian  are  degenerate  and,  strictly  speaking,  the  degenerate  perturbation  theory  rather 
than  the  ordinary  one  should  be  employed.  However,  it  turns  out  that,  though  the  eigenvalues 
(27)-(29)  are  divergent,  the  neutrino  oscillation  probabilities  are  finite  in  the  limits  A  — »  0 
and  A  — >  1.  The  reason  for  this  interesting  behavior  is  a  cancellation  of  divergences  between 
the  eigenvalues  and  the  matrix  elements  of  the  leptonic  mixing  matrix  in  the  calculation  of 
the  evolution  matrix  according  to  Eq.  (10).  In  particular,  in  the  limit  A  — >  0,  the  correct 
vacuum  neutrino  oscillation  probabilities  are  obtained. 

We  shall  now  compare  the  above  results  with  those  existing  in  the  literature.  Equa¬ 
tions  (27)-(32)  have  previously  been  derived  in  Ref.  [35]  and  confirmed  in  Ref.  [36].  Expres¬ 
sions  (33)  and  (34)  are  new.  In  Ref.  [38],  an  expression  for  PA(r  was  found  to  first  order  in 
a,  including  the  0(as  13)  term,  which  can  be  compared  with  the  corresponding  terms  in  our 
Eq.  (34).  We  find  that,  while  our  0(  1)  and  0(a )  terms  coincide  with  those  in  Eq.  (A3)  of 
Ref.  [38],  our  0(as  13)  term  is  quite  different.  In  particular,  we  disagree  with  the  statement 
in  Ref.  [38]  that  to  order  a  the  probability  P/(r  does  not  depend  on  the  CP-violating  phase 
$cp-  The  existence  of  a  term  proportional  to  a  sin  Sep  in  P^r  is  actually  expected,  since  in 
matter  of  constant  density  the  phase  $cp  is  the  s°lc  source  of  T-violation,  and  up  to  the  sign, 
the  T-odd  terms  in  all  three  transition  probabilities  must  be  the  same  [26].  This  is  indeed 
seen  in  Eqs.  (31),  (32),  and  (34).  We  also  note  that  the  term  of  order  asi3  in  Eq.  (A3)  of 
Ref.  [38]  diverges  at  the  atmospheric  resonance  (A  =  1),  while  our  expression  (34)  is  regular 
at  all  physical  values  of  parameters. 

3.2  Vacuum  neutrino  oscillation  probabilities  up  to  second  order  in  ol  and  s13 

The  vacuum  neutrino  oscillation  probabilities  up  to  second  order  in  a  and  s  \ 3  are  given  by 

P™c  =  1  -  c?  sin2  2012  A2  -  4  s23  sin2  A  ,  (35) 

P™c  =  a2  sin2  2012  c\3  A2  +  4  s23  s2Xi  sin2  A  +  2  a  s13  sin  2 0i2  sin  2 023  cos(A  —  <5CP)  A  sin  A  , 

(36) 

P™c  =  a2  sin2  20i2  sl3  A2  +  4  s23  sin2  A  —  2  a  s13  sin  2  012  sin  2  023  cos  (A  —  <5Cp)  A  sin  A  , 

(37) 
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P™c  =  1  -  sin2  2023  sin2  A  +  ac212  sin2  2023  A  sin  2  A 

—  a2  A2  [sin2  2012  c23  +  c22  sin2  2023  (cos  2  A  —  s22)]  +  4  s23  s23  cos  2023  sin2  A 

—  2  a  s33  sin  2012  s23  sin  2023  cos  hep  A  sin  2  A  ,  (38) 

P™c  =  sin2  2023  sin2  A  —  a  c\2  sin2  2 023  A  sin  2  A 

+  a2  sin2  2023  A2  ^ c\2  cos  2  A  —  ^  sin2  2012  sin2  A^ 

—  2  s23  sin2  2023  sin2  A 

+  2  a  si3  sin  2012  sin  2023  (sin  5CP  sin  A  —  cos  2023  cos  hep  cos  A)  A  sin  A  .  (39) 

4  Series  expansion  up  to  first  order  in  a 

In  this  section,  we  expand  the  probabilities  up  to  first  order  in  the  small  parameter  a  while 
keeping  their  exact  dependence  on  013.  One  expects  these  formulas  to  be  useful  for  relatively 
large  values  of  si3.  In  addition,  they  will  correctly  account  for  the  “atmospheric”  resonance 
driven  by  the  parameters  Am2 x  and  0i3. 

4.1  Matter  of  constant  density 

The  eigenvalues  of  the  Hamiltonian  (6)  up  to  first  order  in  a  are 

(40) 

/  \  TD  ^  1  1 

E '2  —  2p 1  2  ^  +  ^  _  2Ci3  a  S 22  ^13  +  —  A  cos  2^13)  ,  (41) 

/\  771  ^  1  1 

E. 3  —  2  ^ 1  2  +  A  + C13)  +  ^—a  —  1  +  Acos20i3)  ,  (42) 

where 

Ci3  =  Y/sin22013  +  {A  —  cos2013)2  .  (43) 

To  calculate  the  probabilities  to  first  order  in  a  we  followed  two  different  approaches:  First, 
we  used  the  Cayley-Hamilton  formalism  as  described  in  Appendix  A. 2,  to  series  expand  the 
evolution  matrix,  and  second,  we  considered  the  constant-density  limit  of  the  expansion  of 
the  evolution  equation  described  in  Appendix  B.l.  Both  methods  gave  the  same  results, 
which  is  a  useful  cross-check  of  our  calculations.  Writing 

P&  =  +  a  P§  +  0(c?) ,  (44) 
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we 


obtain  the  following  expressions  for  the  ve  — >  ve  channel: 
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and  finally,  for  the  vp  — >  uT  channel  we  have 
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+  COS  26*23  cos  $CP 


.  , ,  .  cos  2013  —  A  . 

sm(l  +  A)  A  H - — - sin  C,3A 


C 


13 


X 


sin  Ci  3  A 


(1  +  24A  +  A2)  -  (1  +  A)  sin(l  +  A)  A 


C 


13 


(50) 


In  the  above  formulas,  one  may  identify  the  effective  mixing  angle  in  matter  in  the  1-3  sector 
0'13,  which  is  determined  through  the  expressions 


sin  26*13  = 


sin  2  9 

“4T 


COS  24  = 


cos  2013  —  A 


C 


13 


(51) 


appearing  frequently  in  Eqs.  (45)-(50).  Furthermore,  the  combination  C13A  appearing  as  the 
argument  of  sine  or  cosine  corresponds  to  the  effective  A  in  matter.  In  the  limit  when  013  is 
small,  one  has  C13  ~  A  —  1,  and  expanding  Eqs.  (45)-(50)  up  to  second  order  in  S13  yields 
the  double  expansions  given  in  Sec.  3,  except  for  the  terms  of  order  a2. 

Equations  (40)-(42)  and  (45)-(48)  have  previously  been  derived  in  Ref.  [35]  and  confirmed 
in  Ref.  [36].  Expressions  (49)  and  (50)  are  new. 


4.2  Vacuum  oscillation  probabilities  up  to  first  order  in  a. 


Taking  the  limit  A  — >  0  in  Eqs.  (45)-(50),  it  is  straightforward  to  obtain  the  neutrino  oscil¬ 


lation  probabilities  in  vacuum  to  first  order  in  a: 

Pe(°)vac  =  1  -  sin2  2 013  sin2  A  ,  (52) 

Pe(e1)vac  =  A4  sin2  2di3  sin  2 A  ,  (53) 

PJ4ac=  sm2  20,343  sin2  A,  (54) 

pC0vac  =  _  a4  sin2  2013s23  sin  2A 

+  A  sin  2012Si3C23  sin  2023  (2  sin  5CP  sin2  A  +  cos  $cp  sin  2  A) ,  (55) 

pg))vac  =  c4s  sin2  2  023  Sin2  A  ,  (56) 

4rV “  =  ~Acls  sin2  2023(4  -  44)  sin2A 

+  A  sin  2012  S13C33  sin  2023  (2  sin  5qp  sin2  A  —  cos  5qp  cos  2023  sin  2  A) .  (57) 


5  Series  expansion  up  to  first  order  in  S13 

In  this  section,  we  expand  the  probabilities  up  to  first  order  in  the  small  parameter  s',3  while 
keeping  their  exact  dependence  on  a.  These  formulas  are  expected  to  be  useful  whenever 
neutrino  oscillations  driven  by  the  solar  mass  squared  difference  Am21  are  important.  In 
practice,  this  means  low  neutrino  energies  and  long  baselines. 
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5.1  Matter  of  constant  density 


The  eigenvalues  of  the  Hamiltonian  (6)  np  to  first  order  in  S13  are  given  by 
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(59) 

(60) 


where 
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12 


^cos  20 


2 


(61) 


Note  that  these  eigenvalues  are  independent  of  the  expansion  parameter  s  1 3 ,  which  is  consis¬ 
tent  with  Eqs.  (27)-(29),  where  the  lowest-order  in  S13  corrections  to  the  eigenvalues  appear 
only  at  order  sf3. 

As  in  the  case  of  the  single  expansion  in  a,  we  have  calculated  the  probabilities  using 
two  different  methods:  by  series  expanding  the  evolution  matrix  using  the  Cayley-Hamilton 
formalism  as  described  in  Appendix  A. 2,  and  using  the  constant-density  limit  of  the  single 
expansion  in  .s  1 3  of  the  evolution  equation  described  in  Appendix  B.2.  We  have  found  identical 
results.  Writing  the  probabilities  as 


Poe  =  A?  +  «I3  p§  +  0(4) .  (62) 

we  obtain  the  following  expressions  for  the  ve  — >  ve  channel: 

pee]  =  1  -  Smrl9u  sin2  aCV2A  ,  (63) 

(-/12 

J’i1,=0.  (64) 


The  absence  of  any  first  order  corrections  to  Pee  is  consistent  with  Eq.  (30).  Similarly,  for 
the  ve  — >  vjt  channel  we  find 
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and  for  the  z/M 
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/A-  channel  we  have 
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2  Cl2 


sin2  aC  i2A  —  cos(a:Ci2  +  A  +  a  —  2)  A 


,  cos  26*12  —  —  . 
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sin  2d12  . 

sin  202  3 


x 


1  —  A  —  a  +  Ackc22 

r 
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(cosaC^A  —  cos(A  +  a  —  2)  A)2 
\ 

sinaCi2A  +  sin  (A  +  a  —  2)  A 

sin  aC\2A  +  sin(A  +  a  —  2)  A 


+ 


12 


«AC 


12 


+  sin  Sqp  (1  —  a)  (cosaC^A  —  cos(A  +  a  —  2)  A)  smaC^A  >  . 


(67) 


(68) 


In  this  case,  one  may  identify  in  Eqs.  (63)-(68)  the  effective  “solar”  mixing  angle  in  matter 
$12,  which  is  determined  by 


sin2d,12  = 


sin  2$ 


12 


C 


cos  26*/12  = 


cos2$i2  — 


12 


C 


12 


(69) 


The  combination  aC12A  appearing  as  argument  of  sine  or  cosine  corresponds  to  oscillations 
with  the  “solar”  frequency  in  matter.  Furthermore,  we  note  that  in  the  limit  when  a  is  small, 
one  has  CV2  —  A/ a  —  cos2$i2,  and  expanding  Eqs.  (63)-(68)  up  to  second  order  in  a  yields 
the  double  expansions  given  in  Sec.  3,  except  for  the  terms  of  order  s23. 

Neutrino  oscillations  probabilities  in  matter  of  constant  density  expanded  to  first  order 
in  Si3,  but  exact  in  a ,  presented  in  this  subsection,  have  not  been  previously  published  and 
are  entirely  new. 


5.2  Vacuum  oscillation  probabilities  up  to  first  order  in  S13 

Taking  the  limit  A  — >  0  in  Eqs.  (63) — (68) ,  one  recovers  the  vacuum  probabilities  expanded 
to  first  order  in  S13: 

P"‘rx  =  1  -  sin2  20,2  sin2  aA , 

P«*“  =  0, 

q;;'”"  =  Sill2  2H,  , 4,  Sill2  aA  , 

/') 'c™  —  cos^cp  sin  2012  sit  2023 [sin2  A  —  sin2(l  -  a)  A  +  cos20i2  sin2  aA] 
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(70) 

(71) 

(72) 


(73) 

(74) 


+  -  sin  5cp  sin  26V2  sin  2023 [—  sin  2 A  +  sin  2(1  —  a)  A  +  sin  2a A] , 

Pj^VSjC  =  s212  sin2  29-23  sin2  A  +  c22  sin2  2023[sin2(l  —  ck)  A  —  s22  sin2  aA] , 
pCi)vac  _  CQS  ^cp  sjn  2012  sin  2023  cos  2023[—  sin2  A  +  sin2(l  —  a)  A  +  cos  29 12  sin2  a  A] 

+  -  sin  hCP  sin  2012  sin  2023  [—  sin  2 A  +  sin  2(1  —  a)  A  +  sin  2a A] .  (75) 

6  Qualitative  discussion  and  tests  of  accuracy 

In  this  section,  we  discuss  the  qualitative  behavior  of  the  neutrino  oscillation  probabilities 
and  the  relevance  of  matter  effects  (Sec.  6.1).  We  also  assess  quantitatively  the  accuracy 
of  the  various  expansions  and  compare  each  type  of  expansion  with  an  exact  numerical  cal¬ 
culation  of  the  corresponding  probability  (Secs.  6.2  and  6.3).  We  examine  in  detail  which 
type  of  expansion  is  most  accurate,  depending  on  the  values  of  the  fundamental  parameters 
(especially  a  and  si3),  and  on  the  experimental  configuration  characterized  by  the  neutrino 
energy  E  and  the  baseline  length  L.  In  Sec.  6.4,  we  discuss  issues  related  to  the  applica¬ 
tion  of  our  formulas  to  neutrino  oscillation  experiments.  As  we  have  shown  in  Sec.  2,  the 
probabilities  for  all  oscillation  channels  can  be  obtained  from  the  expressions  for  Pe/Ji  and  P^T 
[see  Eqs.  (16) — (22)] .  Therefore,  we  focus  in  the  following  on  the  ve  — >  v[L  channel;  we  briefly 
comment  also  on  the  other  neutrino  oscillation  channels. 

6.1  The  relevance  of  matter  effects  and  the  probability  Pefl 

Before  considering  the  accuracy  of  the  expansion  formulas,  we  discuss  in  this  subsection 
the  features  and  the  relevance  of  matter  effects  for  the  ve  — >  utJ  oscillation  probability  Pe/l. 
The  discussion  applies  also  to  PeT ;  matter  effects  on  the  probabilities  P/i/t,  PTT,  and  P^T  are 
significantly  weaker  than  those  on  Pefl  and  PeT.  Figure  1  shows  contours  of  Pe/t  for  matter 
of  constant  density  calculated  numerically  from  Eq.  (6)  without  any  approximations,  for  a 
wide  range  of  baseline  lengths  and  neutrino  energies.  Many  features  of  this  figure  can  be 
understood  by  considering  the  expression  for  the  two-flavor  neutrino  oscillation  probability 
in  matter  of  constant  density: 

P  =  sin2  C A2u  ,  C  =  \J sin2  29  +  (cos  29  -  A2u)2  ,  (76) 

where,  in  analogy  to  Eqs.  (25)  and  (26),  we  define  A2l/  =  Am2L/(4E),  A2v  =  2EV/Am2, 
and  9  and  Am2  are  the  generic  two-flavor  neutrino  oscillation  parameters.  From  Eq.  (76) 
it  is  clear  that  for  A2l/  <C  1  one  has  C  ~  1,  and  vacuum  neutrino  oscillations  with  P  ~ 
sin2  29  sin2  A2u  are  recovered.  For  A2u  =  cos  29  one  can  see  in  Eq.  (76)  also  the  MSW 
resonance  [25],  which  leads  to  C  —  sin  29  in  Eq.  (76)  and  to  the  effective  mixing  angle  in 
matter  sin2  29'  =  sin2  2 9/C2  =  1.  The  general  resonance  conditions  for  three  flavors  are  much 
more  complicated  than  those  in  the  two-flavor  case.  However,  due  to  the  hierarchy  of  the 
mass  squared  differences  and  smallness  of  9 is,  one  can  use,  to  a  very  good  approximation,  the 
respective  effective  two-flavor  picture.  This  leads  to  the  “solar”  and  “atmospheric”  resonance 
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Figure  1:  Contours  of  Fe/i  calculated  numerically  without  approximations  for  siir  2ft  13  = 
0.02,  A m31  =  2  •  10_3eV2,  a  =  0.03,  0i2  =  33°,  023  =  45°,  and  5cp  =  0.  Matter  of  constant 
density  p  =  3g/cm3  is  assumed,  and  the  probability  is  averaged  over  a  Gaussian  energy 
resolution  of  1%. 


energies  shown  as  vertical  lines  in  Fig.  1 
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The  “atmospheric”  resonance  is  clearly  visible  in  Fig.  1  at  Eres  ~  8.6  GeV  and  L  >  5500  km. 
The  “solar”  resonance  occurs  at  the  energy  Eres  ~  0.11  GeV.  Note  that  the  solar  resonance  is 
not  as  pronounced  in  Fig.  1  as  the  atmospheric  one,  since  the  neutrino  oscillation  amplitude 
in  vacuum  sin2  2 012  is  already  quite  large. 

Far  enough  to  the  left  of  the  vertical  lines  in  Fig.  1  marking  the  solar  and  atmospheric 
resonance  energies  one  has  vacuum  neutrino  oscillations  with  the  “solar”  parameters  0i2, 
A?7i2i  and  with  the  “atmospheric”  parameters  $13,  A m31,  respectively.  Indeed,  the  typical 
L/ E  pattern  is  clearly  visible  in  Fig.  1  to  the  left  of  the  vertical  lines.  The  diagonal  lines  in¬ 
dicate  the  constant  values  of  L/E  corresponding  to  the  first  solar  and  atmospheric  oscillation 
maxima  in  vacuum,  given  by  the  conditions 
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A  m21L 
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(78) 


respectively.  To  the  right  of  the  vertical  lines  the  probability  is  dominated  by  matter  effects. 
For  A2v  >  1  we  have  C  —  A2v  in  Eq.  (76),  which  has  two  consequences:  First,  since 
v42i/ A2i/  =  VLj 2,  the  oscillation  frequency  becomes  independent  of  energy,  and  second,  the 
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oscillation  amplitude  becomes  suppressed  at  high  energies  because  sin2  2 9  /  A^v  oc  1/E2.  Both 
these  effects  are  apparent  in  Fig.  1  for  neutrino  oscillations  with  the  solar  as  well  as  the 
atmospheric  frequency. 

In  addition  to  the  dependence  on  the  neutrino  energy,  matter  effects  depend  crucially  on 
the  baseline.  From  Eq.  (76)  one  can  see  that  the  two-flavor  oscillation  probabilities  approach 
the  vacuum  ones  for  C A2v  <C  tt/2.  Far  above  the  MSW  resonance  energy,  i.e.  for  A2l/  3>  1, 
one  finds  C  ~  A2v.  Hence,  the  short-baseline  limit  is  C A2v  —  A2v  A2v  =  VL/ 2  <C  7t/2,  leading 
to  vanishing  matter  effects  for  L  <C  Lmat  =  i\ jV  ~  5453  km.  The  horizontal  lines  in  Fig.  1 
indicate  this  baselines  of  the  first  “matter  effect  maximum”  for  large  energies  at  L  =  Lmat, 
and  the  first  “matter  effect  minimum”,  L  =  2 Lmat  ~  10907  km.  At  the  baseline  2Lmat,  where 
A  A  =  VL/ 2  =  7T,  the  terms  proportional  to  a2  and  asi3  in  the  double  expansion  Eq.  (31) 
disappear,  and  only  the  term  proportional  to  s23  survives.  Therefore,  this  baseline  is  especially 
useful  to  measure  S13,  since  ambiguities  due  to  parameter  degeneracies  are  avoided  [38].  A 
dedicated  analysis  of  this  “magic  baseline”  can  be  found  in  Ref.  [46]. 3  Note  that  the  distance 
Lmat  depends  only  on  the  matter  potential  V,  i.e.  the  energy  E  as  well  as  neutrino  masses 
and  mixings  do  not  enter. 

Let  us  now  consider  the  energy  region  well  below  the  MSW  resonance  energy,  where 
A2v  <  1.  The  short-baseline  behavior  of  the  oscillation  probability  in  this  region  can  be 
understood  in  the  two-flavor  picture  by  expanding  Eq.  (76)  for  fixed  A2v  (i.e.,  fixed  L/E) 
assuming  A2v  <C  1  and  A 2uA2v  <C  1.  One  finds 

pmatter  _  pvacuum  ^  ^  Sin2  29  COS  29  sill  A2u  (sin  A2v  —  A2u  COS  A2l)  A2v  .  (79) 

This  relation  shows  that  pmatter  —  pvacuum  decreases  linearly  with  A2v  cx  E,  which  means 
that  the  matter  effects  vanish  for  very  small  energies,  as  mentioned  above.  In  addition,  for 
a  fixed  energy  the  right-hand  side  of  Eq.  (79)  becomes  zero  for  A2v  -C  1  or,  equivalently, 
L  <C  Lvac  =  AnE/Arn2.  Therefore,  the  two  conditions  for  the  smallness  of  matter  effects  due 
to  the  short  baseline  can  be  written  for  all  energies  as  L  -C  Lmin  =  min(Lvac,  Lmat). 

The  relevance  of  matter  effects  is  illustrated  in  Fig.  2.  In  the  upper  left  panel  we  show 
the  probability  difference  |P™tatter  —  -PeV/!'cuum|  f°r  the  case  a  —  0,  such  that  the  oscillations 
due  to  the  “solar”  frequency  are  switched  off  and  only  matter  effects  are  present  which  are 
related  to  the  “atmospheric”  resonance  energy  at  8.6  GeV.  It  can  be  seen  that  the  difference, 
as  well  as  the  relative  difference  (lower  left  panel),  vanish  for  L  <C  Lmin.  The  left  plots  of 
Fig.  2  show  that  the  differences  vanish  also  for  small  energies  and  long  baselines.  This  can 
be  understood  as  follows.  For  E  <C  Eies  the  oscillation  amplitude  in  matter  is  always  close 
to  that  in  vacuum.  For  very  short  baselines,  this  is  also  true  for  the  oscillation  phases.  At 
intermediate  baselines,  the  matter- induced  correction  to  the  oscillation  phase,  though  much 
smaller  than  the  main  “vacuum”  term,  can  become  comparable  to  unity  and  so  cannot  be 
ignored.  However,  at  long  baselines  L  Lvac  the  averaging  regime  sets  in,  which  implies  that 
any  shifts  in  the  oscillation  phase  become  unimportant  and  the  vacuum  oscillations  limit  is 
regained. 

3It  should  be  noted  that  the  quoted  number  for  Lmat  holds  for  a  constant  matter  density  of  p  =  3g/cm3. 
It  will  differ  for  larger  matter  densities  like  those  in  the  mantle  of  the  Earth,  or  even  more  drastically 
in  astrophysical  or  cosmological  applications  of  neutrino  oscillations.  In  addition,  the  approximation  of 
constant  matter  density  is  then  frequently  not  justified.  For  a  realistic  Earth  matter  density  profile  one  finds 
2Tmat  ~  7250  km  [46]. 
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Figure  2:  Contours  of  |P“atter  -  P”m|  (upper  panels),  and  |P”atter  -  Pe”m|/P“ “tter 
(lower  panels),  where  p™atter  (PJacuum)  is  the  exact  oscillation  probability  in  matter  (vacuum) 
for  sin2  26^  =  0.05,  Am2,  =  2  ■  10-3  eV2,  and  023  =  45°.  In  the  left  panels  a  =  0,  whereas  in 
the  right  panels  a  =  0.03,  0\2  =  33°,  and  <5cp  =  0.  Matter  of  constant  density  p  =  3g/cm3 
is  assumed,  and  the  probabilities  are  averaged  over  a  Gaussian  energy  resolution  of  2%.  The 
straight  lines  indicate  the  relevant  resonance  energies,  first  oscillation  maxima,  and  the  first 
matter  effect  maximum  (see  also  Fig.  1). 


The  above  discussion  describes  the  regions  where  |Pe™atter  —  Pevacuum  vanishes  or  becomes 
small  for  a  =  0  in  the  left  plots  of  Fig.  2.  One  can  also  understand  in  this  way  why  the 
difference  |Pe“atter  —  P^cuum|  disappears  more  slowly  towards  low  E  and  L  along  the  lines 
of  constant  P/P  in  the  upper  left  plot  of  Fig.  2.  The  structures  seen  in  the  plot  emerge 
from  the  matter  dependent  phase  shifts  of  oscillation  probabilities  with  more  or  less  equal 
amplitudes.  This  behavior  can  be  understood  in  the  two- flavor  picture  from  Eq.  (79),  which 
implies  that  pmatter  —  pvacuum  grows  linearly  with  P  along  the  lines  of  constant  /S.2V  or, 
equivalently,  constant  P/P.  This  growth  is  modulated  in  the  A2„-direction  by  the  factor 
sin  A2u  (sin  A2i,  —  A2y  cos  A2i/),  which  describes  nicely  the  details  of  the  structures  extending 
along  lines  of  constant  A2u  towards  lower  energies. 

The  right  plots  of  Fig.  2  show  the  probability  difference  |Pe™atter  —  -P™cuum|  and  the 
relative  difference  for  the  case  a  =  0.03.  This  figure  contains  the  structures  stemming  from 
the  atmospheric  resonance  energy,  which  were  already  present  in  the  a  —  0  case  (left  plots). 
In  addition,  similar  structures  show  up  around  the  solar  resonance  energy,  reaching  to  lower 
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Figure  3:  Rows  1,  2,  and  3:  absolute  errors  of  the  three  types  of  expansions  for  PefJ  as 
functions  of  neutrino  energy  E  and  baseline  length  L,  for  Am2,  =  2  ■  10~3eV2,  a  =  0.03, 
#12  =  33°,  #23  =  45°,  <5cp  =  0,  and  several  values  of  sin2  2#13.  Matter  of  constant  density 
p  =  3g/cm3  is  assumed.  Row  4  shows  contours  of  Pe/t  calculated  numerically. 


energies  and,  due  to  the  larger  effective  mixing  angle,  also  to  shorter  baselines.  The  right 
plots  of  Fig.  2  show  that  in  the  realistic  case  of  two  mass  squared  differences  matter  effects 
are  rather  important  and  have  to  be  taken  into  account  in  a  large  domain  of  the  physically 
interesting  parameter  space. 

6.2  Comparing  the  accuracy  of  the  three  types  of  expansions 

In  order  to  test  the  accuracy  of  our  analytic  expressions,  we  shall  now  compare  the  values 
of  Pefl  obtained  from  the  expansion  formulas  Eq.  (31)  for  the  double  expansion,  Eqs.  (47) 
and  (48)  for  the  single  expansion  in  a,  and  Eqs.  (65)  and  (66)  for  the  single  expansion  in 
Si3  with  the  exact  values  of  Pefl  calculated  numerically.  Figure  3  shows  the  absolute  errors 
|pexpanslon  _  pex act  | ,  anc[  p^g  4^  the  relative  errors  P^xpansion  _  pexacty pexact  for  the  three 

types  of  expansions  as  functions  of  the  neutrino  energy  E  and  baseline  length  L  for  various 
values  of  sin22#13.  For  reference,  we  display  in  the  lowest  row  of  graphs  in  each  of  these 
figures  also  the  probability  PC)1  itself,  calculated  numerically  without  any  approximations. 

First,  we  observe  that  in  general  the  absolute  errors  shown  in  Fig.  3  are  rather  small — 
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Figure  4:  Rows  1,  2,  and  3:  relative  errors  of  the  three  types  of  expansions  for  Pe[l  as 
functions  of  neutrino  energy  E  and  baseline  length  L,  for  Am2,  =  2  •  10~3eV2,  a  =  0.03, 
d\2  =  33°,  d23  =  45°,  <5cp  =  0,  and  several  values  of  sin2  2 d13.  Matter  of  constant  density 
p  =  3g/cm3  is  assumed.  Row  4  shows  contours  of  Pe M  calculated  numerically. 

at  the  0.1%  level  in  a  large  part  of  the  parameter  space,  whereas  the  relative  errors  shown 
in  Fig.  4  are  considerably  larger,  due  to  the  smallness  of  the  probability  itself.  Next,  we 
note  that  the  series  expansions  in  a,  be.,  the  double  expansion  in  a  and  S13  and  the  single 
expansion  in  a,  are  only  valid  for 

A?n2  L  L 

aA  =  —  21  -Cl,  or  —  -C  104  km/GeV ,  (80) 

4  E  E 

be.,  far  below  the  first  solar  maximum.  The  obvious  reason  is  that  expanding  terms  of  the 
type  sinaA  is  only  valid  if  aA  is  small,  and  hence,  these  two  types  of  expansion  cannot 
account  for  neutrino  oscillations  with  the  “solar”  frequency.  Note  that  for  S13  =  0  the  single 
expansion  in  a  gives  Pefl  =  0,  since  in  that  case  the  lowest-order  term  of  Pe ^  is  proportional 
to  a2,  and  our  single  expansion  in  a  only  contains  terms  up  to  first  order  in  a.  This  explains 
why  at  very  small  values  of  S13  the  double  expansion  (which  includes  terms  of  second  order 
in  a)  is  more  accurate  than  the  single  expansion  in  a.  In  contrast,  the  neutrino  oscillations 
with  the  “solar”  frequency  (for  which  the  condition  in  Eq.  (80)  is  violated)  are  in  general 
very  well  accounted  for  by  the  single  expansion  in  S13,  since  it  retains  the  exact  dependence 
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of  the  probability  on  a.  This  expansion  is  the  best  one  for  relatively  small  values  of  .s  1 3  and 
large  values  of  L/E. 

Figures  3  and  4  allow  us  to  put  the  above  observations  on  a  more  quantitative  basis.  In 
the  region  of  L/E  given  in  Eq.  (80)  (lower-right  parts  of  the  graphs),  where  the  oscillations 
are  mainly  driven  by  the  “atmospheric”  mass  squared  difference  the  double  expansion 

and  the  single  expansion  in  a  work  rather  well.  For  not  too  large  values  of  sin2  2 013  the 
double  expansion  is  most  accurate,  with  absolute  errors  at  the  0.1%  level  and  relative  errors 
not  exceeding  1%  (5%)  for  sin2  26h3  =  0.001  (0.01).  However,  for  values  of  sin2  26h3  close  to 
the  current  upper  bound  the  single  expansion  in  a  becomes  better,  with  relative  errors  smaller 
than  1%.  The  reason  for  this  is  that  for  large  values  of  sin2  26,13  the  neutrino  oscillations  driven 
by  the  “atmospheric”  frequency  and  mixing  angle  6 13  completely  dominate  the  probability, 
and,  in  addition,  for  relatively  large  values  of  L  and  energies  close  to  10  GeV,  the  atmospheric 
resonance  becomes  important.  Since  the  single  expansion  in  a  is  exact  in  .s  1 3 ,  it  describes  the 
case  of  relatively  large  s33  well,  and  the  atmospheric  resonance  is  also  correctly  accounted 
for.  At  the  same  time,  the  accuracy  of  the  double  expansion  becomes  slightly  worse  near 
the  first  atmospheric  maximum  and  the  resonance.  As  was  pointed  out  in  Ref.  [36],  for  large 
values  of  s33  the  accuracy  of  the  double  expansion  can  be  improved  by  replacing  the  term 
proportional  to  s23  in  Eq.  (31)  by  the  term  given  in  Eq.  (47)  as  the  zeroth  order  term  in  the 
single  expansion  in  a.  For  a  =  0  this  term  describes  the  probability  Peil  exactly  to  all  orders 
in  si3. 

As  can  be  seen  from  Figs.  3  and  4,  the  single  expansion  in  Si3  gives  a  rather  poor  de¬ 
scription  of  the  region  of  L/E  defined  in  Eq.  (80).  The  reason  is  that  with  only  terms  of  first 
order  in  s33  it  is  not  possible  to  obtain  a  correct  description  of  neutrino  oscillations  driven 
by  A m31  and  On.  This  is  also  reflected  by  the  fact  that  the  lowest-order  in  s33  terms  in 
the  eigenvalues  of  the  Hamiltonian  are  (R(s23)  [see  Eqs.  (27)-(29)].  On  the  other  hand,  the 
single  expansion  in  s33  is  rather  accurate  for  large  values  oi  L/E  violating  the  condition  (80) 
(upper-left  parts  of  the  graphs)  and  relatively  small  values  of  6 13:  for  sin2  2 d13  <  0.001  the 
accuracy  is  typically  better  than  1%. 

To  conclude  this  subsection,  we  show  in  Fig.  5  which  type  of  expansion  provides  the  most 
accurate  expression  for  PeM,  depending  on  the  values  of  the  expansion  parameters  a  and  Si3 
and  for  a  number  of  fixed  values  of  neutrino  energy  E  and  baseline  length  L.  These  plots 
change  very  little  when  the  fundamental  neutrino  parameters  A m31,  012,  023,  and  <5cp  are 
varied  within  their  allowed  ranges,  and  also  when  one  switches  over  to  the  other  neutrino 
oscillation  channels.  As  expected,  one  observes  as  a  general  trend  that  the  single  expansion 
in  a  is  best  for  small  a  and  large  Si3,  whereas  the  single  expansion  in  Si3  is  best  for  small 
Si3  and  large  a.  The  double  expansion  is  most  accurate  in  a  region  where  a  and  s33  are  of 
comparable  order  of  magnitude. 

In  agreement  with  the  discussion  related  to  Figs.  3  and  4,  we  find  from  Fig.  5  that  in  the 
low-energy  regime  P  ~  0.1  GeV  the  single  expansion  in  Si3  is  most  accurate  in  almost  the 
entire  a-0 13  plane.  For  values  oiL/E  satisfying  the  condition  (80)  and  energies  larger  than  a 
few  GeV  the  double  expansion  is  most  accurate  in  a  large  fraction  of  the  a-0 33  plane;  only  for 
values  of  s33  close  to  the  current  upper  bound  does  the  single  expansion  in  a  become  better. 
Furthermore,  we  note  that  for  E  ~  10  GeV  and  L  >  5500  km  the  atmospheric  resonance  is 
important,  which  leads  to  a  better  accuracy  of  the  single  expansion  in  a  (see  the  rightmost 
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Figure  5:  The  plot  shows  which  type  of  expansion  for  Pe/J  is  most  accurate,  depending  on 
a  and  S13  and  for  different  values  of  neutrino  energy  E  and  baseline  length  L.  The  values 
of  the  fundamental  neutrino  parameters  used  are  Am^  =  2  •  10~3eV2,  6*12  =  33°,  $23  =  45°, 
and  <5cp  =  0.  Matter  of  constant  density  p  =  3g/cm3  is  assumed.  The  dotted  lines  indicate 
the  3cr  allowed  range  of  the  parameter  a. 


panel  in  the  third  row  of  Fig.  5). 

We  conclude  that  the  double  expansion  in  both  a  and  S13  is  most  accurate  in  a  wide 
region  of  the  parameter  space,  where  a  and  S13  are  roughly  of  the  same  order  of  magnitude. 
The  single  expansion  in  S13  has  to  be  used  whenever  neutrino  oscillations  with  the  solar 
frequency  are  important  (e.g.,  in  the  low-energy  regime),  whereas  the  single  expansion  in  a 
is  most  accurate  for  values  of  s  !3  close  to  the  upper  bound,  or  in  cases  where  the  atmospheric 
resonance  is  important. 

6.3  Accuracy  of  the  double  expansion 

Motivated  by  the  fact  that  the  double  expansion  is  most  accurate  in  a  wider  region  of  the 
parameter  space  than  the  single  expansions,  we  present  in  this  subsection  some  more  accuracy 
tests  for  it.  Figures  6  and  7  show  the  relative  errors  for  Pefl  given  in  Eq.  (31)  and  P/lT  given  in 
Eq.  (34),  respectively,  for  neutrinos,  antineutrinos,  and  normal  and  inverted  hierarchies.  We 
note  that  the  relative  accuracy  of  PeT  and  1  —  Pee  is  very  similar  to  the  one  of  Pefl  (Fig.  6), 
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Figure  6:  Relative  errors  of  the  double  expansion  for  Pefl  for  neutrinos  (rows  1  and  2),  and 
PfTjx  for  antineutrinos  (rows  3  and  4),  for  normal  hierarchy  (NH)  and  inverted  hierarchy  (IH). 
The  values  of  the  fundamental  neutrino  parameters  used  are  =  2-  10_3eV2,  a  =  0.03, 

6 12  =  33°,  6*23  =  45°,  hep  —  0,  and  several  values  of  sin22#13.  Matter  of  constant  density 
p  =  3  g/ cm3  is  assumed. 

and  the  accuracy  of  and  PTT  is  similar  to  that  of  P^T  (Fig.  7). 

In  the  region  of  L/E  defined  in  Eq.  (80),  the  accuracy  of  Pefl  is  roughly  between  1%  and 
5%.  It  should  be  noted  that  the  absolute  errors  of  the  double  expansion  of  Pefl  depend  very 
weakly  on  the  choice  of  the  fundamental  neutrino  parameters  within  their  allowed  ranges,  and 
they  are  always  very  small,  similar  to  the  Erst  row  of  plots  in  Fig.  3.  However,  the  relative 
errors  shown  in  Fig.  6  are  rather  sensitive  to  the  fundamental  parameters  values  in  the  region 
of  low  L/E.  The  reason  is  that  the  probability  Pefl  itself  is  typically  very  small  in  this  region, 
being  as  tiny  as  10~7,  or  even  smaller.  Hence,  the  relative  error,  being  a  ratio  of  two  small 
numbers,  is  quite  sensitive  to  variations  of  the  parameters  (see,  e.g.,  the  second  column  in 
Fig.  6).  In  contrast,  P^T  is  rather  large,  because  of  the  oscillations  in  “atmospheric”  channel 
driven  by  the  large  Am2 x  and  maximal  or  nearly  maximal  mixing  sin2  2023  —  1.  This  is  also 
clear  from  Eq.  (34),  since  PMT  has  a  term  of  zeroth  order  in  both  a  and  Si3,  while  Pefl  from 
Eq.  (31)  is  of  second  order  in  these  parameters.  Hence,  in  the  case  of  PllT ,  the  relative  errors 
are  of  the  same  order  as  the  absolute  errors,  which  are  always  very  small.  The  errors  shown 
in  Fig.  7  are  smaller  than  0.1%  for  sin2  2d13  <  0.01  and  smaller  than  1%  for  sin2  2d13  <  0.1. 
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Figure  7:  Relative  errors  of  the  double  expansion  for  P^T  for  neutrinos  (rows  1  and  2),  and 
Pjj: f  for  antineutrinos  (rows  3  and  4),  for  normal  hierarchy  (NH)  and  inverted  hierarchy  (IH). 
The  values  of  the  fundamental  neutrino  parameters  used  are  =  2-  10_3eV2,  a  =  0.03, 

9 12  =  33°,  6*23  =  45°,  5cp  =  0,  and  several  values  of  sin22#13.  Matter  of  constant  density 
p  =  3  g/ cm3  is  assumed. 


6.4  Probability  expansions  and  relation  to  experiments 

The  discussion  of  the  accuracy  of  the  different  expansions  can  be  used  to  identify  the  best 
set  of  equations  for  a  given  neutrino  oscillation  experiment.  The  first  question  is,  however, 
whether  matter  effects  are  relevant  or  if  the  much  simpler  vacuum  probabilities  can  be  used. 
The  absolute  and  relative  differences  between  matter  and  vacuum  probabilities  vanish,  as 
discussed  in  Sec.  6.1,  for  sufficiently  low  energies  or  for  short  enough  baselines.  The  discussion 
showed  that  in  experimentally  interesting  cases  matter  effects  are  almost  always  relevant  for 
the  ranges  of  L  and  E  that  we  considered.  If  L/E  is  such  that  oscillations  can  be  observed, 
then  for  fixed  L/E  the  matter  effects  are  negligible  only  at  very  small  E  (and  consequently 
small  L).  An  example  where  matter  effects  can  be  ignored  at  the  percent  level  is  given  by 
the  reactor  ve  disappearance  experiments  with  energies  of  a  few  MeV  and  baselines  up  to  a 
few  kilometers. 

In  order  to  identify  the  best  expansion  for  a  given  experiment,  one  must  take  into  account 
that  the  analyses  of  real  experiments  are  based  on  event  rates  and  not  on  probabilities.  It 
is  useful  to  distinguish  here  between  “point  sources”  (accelerator  beams,  reactors,  the  Sun, 
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supernovae,  etc.)  where  the  observed  neutrinos  come  from  one  point-like  region  in  space,  and 
“extended  sources”  ( e.g .,  the  atmosphere)  where  the  observed  neutrinos  come  from  many 
different  positions  in  space.  An  event  rate  based  analysis  of  point  sources  implies  essentially, 
for  a  given  energy,  an  additional  factor  1  / L2  in  order  to  account  for  the  reduction  of  the 
flux  as  a  function  of  the  baseline  L.  In  addition,  the  detection  cross  section  is  typically 
proportional  to  Er ,  where  1  <  r  <  2,  depending  on  the  energy  range  and  on  the  detection 
process.  In  the  two-flavor  approximation,  when  going  from  the  first  oscillation  maximum  at 
/^2u  =  7r/2  to  the  nth  oscillation  maximum  at  A2l/  =  (n  —  1/2)7t  this  implies  that  the  event 
rates  drop  by  a  factor  l/(2n  —  l)2  if  the  baseline  is  increased,  or  by  a  factor  1/(2 n  —  l)r  if  the 
energy  is  decreased.  Hence,  e.g.,  for  a  fixed  energy  the  rate  in  the  second  (third)  maximum 
is  already  reduced  by  a  factor  1/9  (1/25).  This  explains  why  most  current  proposals  for  long 
baseline  oscillation  experiments  aim  at  the  first  oscillation  maximum  of  the  “atmospheric” 
oscillations.  The  oscillations  driven  by  the  “solar”  frequency  are  then  a  sub-leading  effect, 
and  in  most  cases  the  double  expansion  works  very  well.  The  overall  relative  precision  of  the 
double  expansion  is  then  typically  in  the  range  of  a  percent,  or  at  worst,  up  to  a  few  percent 
in  this  case.  However,  if  $13  is  close  to  its  current  upper  limit,  the  single  expansion  in  a  is 
even  more  precise  than  the  double  expansion.  On  the  other  hand,  if  #13  is  tiny,  well  below 
the  sensitivity  limits  of  any  planned  accelerator  experiment,  the  single  expansion  in  S13  is 
more  precise  than  the  double  expansion,  especially  for  low  energies.  Note,  however,  that  the 
double  expansion  has  often  a  sufficient  precision  even  when  it  is  not  the  best  expansion. 

The  discussion  of  point  sources  does  not  change  much  for  experiments  which  aim  at 
higher  oscillation  maxima.  The  1/(2 n  —  l)2  flux  factor  due  to  the  beam  divergence  and/or 
the  l/(2n  —  l)r  factor  due  to  the  energy  dependence  of  the  detection  cross  sections  reduce 
the  observed  event  rates  and  limit  all  proposals  to  the  first  few  oscillation  maxima  at  best. 
The  double  expansion  still  works  quite  well  for  such  proposals,  as  long  as  the  sub-leading 
oscillations  governed  by  the  “solar”  frequency  stay  in  the  linear  regime,  where  an  expansion 
in  a  makes  sense.  This  leads  to  the  condition  2n  —  1  <C  or1  ~  40,  which  is  fulfilled  for  the  first 
few  maxima.  Higher  values  of  n  are  currently  not  proposed,  since  the  flux  and/or  cross  section 
would  drop  by  a  large  factor,  which  would  make  the  detectors  and  sources  unaffordable.  For 
largest  (smallest)  $13  the  single  expansion  in  a  (S13)  works  again  numerically  even  better  then 
the  double  expansion,  just  like  in  the  case  of  the  first  oscillation  maximum. 

The  discussion  is  somewhat  different  for  sources  which  are  not  point-like.  In  that  case, 
there  is  in  general  no  1  /L2  suppression  of  the  flux,  and  the  observed  event  rates  are  not  dom¬ 
inated  by  the  first  oscillation  maximum.  An  important  example  are  atmospheric  neutrinos, 
for  which  a  wide  range  of  baselines  and  energies  contributes:  the  energy  window  of  practical 
interest  ranges  from  a  100  MeV  to  a  few  tens  of  GeV,  while  the  baseline  lies  between  10  km 
and  104  km.  This  wide  range  of  energies  and  baselines  makes  clear  that  none  of  the  discussed 
expansions  covers  the  full  parameter  range  at  the  percent  level.  In  general  the  double  ex¬ 
pansion  works  quite  well  for  small  or  moderate  L/E ,  while  the  S13  expansion  tends  to  work 
best  for  larger  L/E  values.  At  the  same  time,  close  to  the  atmospheric  resonance  or  if  the 
mixing  angle  $13  is  near  its  current  upper  bound  and  for  not  too  large  L/E ,  the  a  expansion 
is  expected  to  be  the  most  relevant  one. 

In  general,  for  all  types  of  experiments  and  for  any  given  values  of  L,  E,  and  #13  the  most 
accurate  expansion  can  readily  be  found  with  the  help  of  Fig.  5. 
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7  Summary  and  conclusions 


In  this  paper  we  have  presented  three  different  sets  of  approximate  formulas  for  the  prob¬ 
abilities  of  neutrino  oscillations  in  matter  and  in  vacuum.  We  have  shown  that  in  general 
the  probabilities  for  all  possible  oscillation  channels  for  neutrinos  as  well  as  for  antineutrinos 
can  be  obtained  from  just  two  independent  probabilities  by  using  unitarity  and  a  symmetry 
related  to  the  “atmospheric”  mixing  angle  623  [see  Eq.  (14)].  One  possible  choice  for  these 
two  probabilities  is  Pe/l  and  P^T.  We  have  derived  the  expressions  for  the  neutrino  oscilla¬ 
tion  probabilities  in  matter  of  constant  density,  expanded  in  terms  of  the  small  parameters 
a  =  S13  or  in  both  of  them.  Below  we  summarize  the  main  features  of  these 

expansions: 

•  Double  expansion  up  to  second  order  in  a  and  s13  [Sec.  3.1,  Eqs.  (30) — (34)] 

The  neutrino  oscillation  probabilities  are  expanded  in  both  small  parameters  up  to 
second  order.  In  general,  these  expressions  are  valid  for  aA  =  Am^L/^AE)  <C  1 
or  L/E  <C  104km/GeV,  i.e.,  if  neutrino  oscillations  due  to  the  solar  mass  squared 
difference  are  not  important.  The  accuracy  of  the  approximation  is  good  for  a 

wide  range  of  the  parameters.  Typically,  the  relative  errors  of  Pejl ,  PeT,  and  1  —  Pee  are 
between  1%  and  5%,  the  relative  errors  of  P/t/t ,  P/tT ,  and  PTT  and  the  absolute  errors  for 
all  probabilities  are  of  the  order  0.1%. 

•  Single  expansion  up  to  first  order  in  ol  [Sec.  4.1,  Eqs.  (45)-(50)] 

The  neutrino  oscillation  probabilities  are  expanded  in  a,  but  the  exact  dependence  on 
S13  is  retained.  Like  in  the  case  of  the  double  expansion,  these  formulas  are  valid  in 
the  region  where  the  oscillations  driven  by  the  solar  mass  squared  difference  are  not 
important:  aA  =  Am^L/  (AE)  <  1,  or  L/L  <C  104km/GeV.  The  accuracy  is  better 
than  the  one  of  the  double  expansion  for  values  of  S13  >  0.1  close  the  current  upper 
bound,  or  if  the  atmospheric  resonance  is  important.  For  instance,  for  a  matter  density 
p  =  3  g/cm3  and  A m31  =  2  •  10-3  eV2  this  is  the  case  for  E  ~  10  GeV  and  L  >  5500  km. 

•  Single  expansion  up  to  first  order  in  s13  [Sec.  5.1,  Eqs.  (63) — (68)] 

The  neutrino  oscillation  probabilities  are  expanded  in  s  1 3 ,  but  the  exact  dependence 
on  a — and  therefore  on  Am^ — is  retained.  Hence,  these  expressions  are  useful  in  the 
region  where  the  oscillations  due  to  the  solar  mass  squared  difference  are  relevant: 
aA  =  Am|1L/(4E)  >  1,  or  L/E  >  104km/GeV.  In  particular,  this  is  the  case  for  low 
energies  E  ~  0.1  GeV  and  L  >  103  km.  For  values  oi  L/E  outside  the  “solar”  regime 
this  type  of  expansion  is  only  useful  for  very  small  values  of  S13  <  10-4  -j- ICG3. 

•  Vacuum  limit  of  the  expansions  [Secs.  3.2,  4.2,  5.2] 

The  vacuum  limit  for  each  type  of  expansion  is  valid  if,  in  addition  to  the  requirements 
ensuring  the  validity  of  a  given  type  of  expansion,  matter  effects  can  be  neglected  (see 
discussion  in  Sec.  6.1).  This  usually  implies  low  energies  or  short  baselines.  We  find 
that  in  many  realistic  cases  matter  effects  are  of  the  order  of  a  few  percent.  Hence, 
if  an  accuracy  at  that  level  is  required,  one  should  employ  the  formulas  for  neutrino 
oscillations  in  matter. 
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In  conclusion,  we  have  presented  a  collection  of  formulas  for  three-flavor  neutrino  oscilla¬ 
tion  probabilities  by  deriving  expansions  in  small  parameters.  We  have  performed  a  detailed 
analysis  of  the  accuracy  of  these  expansions  and  determined  the  parameter  regions  where 
they  are  most  accurate.  The  expansions  of  the  neutrino  oscillation  probabilities  in  matter 
of  constant  density  are  useful  for  the  analytical  understanding  of  the  physics  of  future  neu¬ 
trino  oscillation  experiments.  Furthermore,  we  have  also  presented  expansion  formulas  for 
the  neutrino  oscillation  probabilities  in  arbitrary  matter  density  profiles  (see  Appendix  B), 
which  can  be  applied  to  a  large  class  of  problems. 
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A  Hamiltonian  diagonalization  approach 

In  this  appendix,  we  present  the  details  of  the  approximate  diagonalization  of  the  effective 
Hamiltonian  of  the  neutrino  system  in  the  case  of  matter  of  constant  density. 

A.l  Diagonalizing  the  Hamiltonian  perturbatively 

In  order  to  derive  the  double  expansions  given  in  Sec.  3,  we  write  the  Hamiltonian  of  Eq.  (6) 
as 

H  ~^ko23UsMUlOls,  (Al) 

where  M  =  Oi3Oi2  diag  (0,  a,  1)  Oj2Of3  +  diag(A,0,0),  and  the  matrices  Oij  and  Us  have 
been  defined  after  Eq.  (1).  The  matrix  M  can  be  explicitly  written  as  M  =  (Am^l/2E)~1  H' 
by  setting  <5cp  =  0  in  H ',  which  is  given  in  Eq.  (Bl)  below.  First,  we  diagonalize  the  matrix 
M  by  M  =  WMW'  with  M  =  diag(A1,  A2,  A3)  and  W  being  a  unitary  diagonalizing  matrix. 
This  diagonalization  is  performed  by  using  perturbation  theory  up  to  second  order  in  the 
small  parameters  a  and  si3,  i.e.,  we  write  M  =  +  ibC1)  +  M^2\  where  M ^  (M^) 

contains  all  terms  of  first  (second)  order  in  a  and  si3.  One  finds 

M(0)  =  diag  (A,  0, 1)  =  diag(Ai°},  A^0),  A^0)) ,  (A2) 
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M(1) 


M(2) 


asf.2 

OtSi2C\2  S13  \ 

atsi2ci2 

«C?2  0 

*13 

0  0  / 

s2 

*13 

0  — a<s13*i2 

0 

0  —  CCS13S12C12 

— CCS13SJ2 

— CCS13S12C12  — *13 

(A3) 

(A4) 


For  the  eigenvectors  we  write  v*  =  vf  '1  +  vp\  and  since  is  diagonal  at  zeroth  order, 

we  have  =  e*.  Then,  the  first  and  second  order  corrections  to  the  eigenvalues  are  given 
by 


A,(1) 

Af 


M 


(i) 


M 


(2) 


£ 


M-P 


x (°)  x(0)  ’ 

jfr  \  Aj 


and  the  corrections  to  the  eigenvectors  are  calculated  by 


= 


(2) 

V)  = 


£ 


M 


(i) 


v 


A(0)  _  A(0)  eJ  ’ 


jW  *  '\7 

1 


£ 


x(°)  \(0) 

j+i  \  Aj 


M[?  +  (M^v?)  -  A,(1)  («<‘>) 


e3- 


(AS) 

(A6) 

(A7) 

(AS) 


The  mixing  matrix  in  matter  is  given  by  U'  =  O23U5W  with  W  =  {yi,V2,v^))  and  the 
eigenvalues  of  the  Hamiltonian  [see  Eqs.  (27)-(29)]  are  obtained  as  =  [Am^/ (2E)]\i 
(*  =  1,2,3). 

In  our  paper,  we  do  not  in  general  order  the  eigenvalues  according  to  their  magnitude. 
Such  an  ordering  would  create  problems  as  one  would  have  to  re-label  the  eigenvalues  upon 
passing  through  each  of  the  two  MSW  resonances.  The  ordering  is  actually  unimportant  if 
one  is  careful  to  assign  the  correct  eigenvector  to  each  eigenvalue. 

We  also  reiterate  the  point  discussed  at  the  end  of  Sec.  3.1:  despite  the  fact  that  in  the 
case  of  the  double  expansion  the  eigenvalues  as  well  as  certain  entries  of  the  leptonic  mixing 
matrix  in  matter  are  divergent  at  A  —a  0  and  A  —a  1,  the  neutrino  oscillation  probabilities 
are  finite  in  these  limits.  In  particular,  the  correct  vacuum  probabilities  are  recovered  in  the 
limit  A  —a  0.  The  mentioned  divergences  are  of  no  concern  to  us,  since  we  are  interested 
in  oscillation  probabilities  rather  than  in  eigenvalues  or  the  mapping  between  the  mixing  in 
matter  and  in  vacuum,  as  was  the  case,  e.g.,  in  Ref.  [36]. 


A. 2  The  Cayley-Hamilton  formalism 


In  order  to  derive  a  series  expansion  for  the  evolution  matrix  S  in  small  parameters,  we  can 
use  the  following  formula  from  Refs.  [12,14]: 

S(t,t0)  =  e~lHL 

g-i-Ei  L 

= -  \E2E3I3  —  ( E2  +  Ez)H  +  H 2] 

(E1  -E2)(E1  -E3)  1  3  3  1  2  3)  J 

g— 1E2L 

+  7 - — — - lE^ls  -  (Ex  +  E3)H  +  H 2] 

(E2-El)(E2-E3)V  13 )  ^  J 

p—iE^L 

+  — - — — - —  \E!E213  -  (E1  +  E,)H  +  H2]  ,  (A9) 

(E3-E1)(E3-E2)  1  J  v  J  v  ; 

where  El  is  the  Hamiltonian  in  the  flavor  basis  as  given  in  Eq.  (6)  and  Ei  (i  =  1,2,3)  are 
the  eigenvalues  of  this  Hamiltonian.  Inserting  series  expansions  of  the  eigenvalues  [see,  e.g., 
Eqs.  (40)-(42)  and  (58)-(60)]  and  the  Hamiltonian  E[  into  Eq.  (A9),  it  is  straightforward 
to  obtain  an  expression  for  Sag,  and  subsequently,  for  the  neutrino  oscillation  probabilities 
PaP=\Spa\2. 


B  Perturbative  expansion  of  the  evolution  equation  for  arbitrary 
matter  density  profiles 


In  this  appendix,  we  present  the  details  of  the  perturbative  expansion  of  the  neutrino  evolution 
equation  in  the  case  of  matter  with  an  arbitrary  density  profile.  We  also  give  the  formulas 
relevant  for  the  particular  case  of  matter  of  constant  density. 

Consider  the  neutrino  evolution  matrix  S(t,to)  defined  in  Eq.  (9).  It  is  convenient  to 
rotate  it  by  the  angle  623  according  to  Eq.  (13).  The  evolution  matrix  in  the  rotated  ba¬ 
sis  S'(t,  t0)  satisfies  the  evolution  equation  i(d/dt)S'(t,  to)  =  H'{t)S'{t}t0)  with  the  initial 
condition  S' (to,  to)  =  t  and  the  Hamiltonian  H'(t)  given  by  [28] 


H\t) 


Here 


s12c13^21  +  s13^31  +  V  (t) 

'S12C12C13A21 

Si3Ci3  (A31  —  s^2A2i)  el5cp 


'S12C12C13A21  S13C13  (A31  —  <s^2 A21 )  e  lScp 

Ci2A2l  — Si2Ci2Si3e  iScpA2i 

— 'Si2Ci2Si3el5cpA2i  c13^31  +  S12S13^21 


=  Am21 
2  E 


A31 


2  E 


It  is  important  to  notice  that  H'(t)  [and  hence,  S'(t,  to)]  does  not  depend  on  the  mixing  angle 
$23-  This  is  a  consequence  of  the  specific  parameterization  of  the  leptonic  mixing  matrix  U 
in  Eq.  (1),  in  which  the  matrix  O23  is  the  leftmost  one,  and  the  fact  that  the  matrix  of  the 
matter- induced  potentials  in  the  effective  Hamiltonian  H(t)  commutes  with  O 23.  Once  the 
matrix  S'(t,t0)  is  found,  one  can  use  Eq.  (13)  to  rotate  back  to  the  original  flavor  basis. 
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It  is  convenient  to  decompose  the  Hamiltonian  H'(t )  as  H'(t )  =  H^t)  +  H),  where  H^t) 
is  of  zeroth  order  in  some  expansion  parameter  and  H)  is  the  remainder.  Accordingly,  the 
evolution  matrix  can  be  written  as  S'(t,t0)  =  S'0(t,t0)S[(t,t0),  where  S0(t,t0 )  satisfies 

t0)  =  Ho(t)S'0(t,  t0) ,  Sq  (t0,  t0)  =  l.  (B3) 

Then  the  matrix  S[(t,t0)  satisfies 

i^(Mo)  =  Si  (to,  to)  =  1-  (B4) 

Up  to  now,  no  approximations  have  been  made.  Next,  we  shall  find  the  evolution  matrix 
perturbatively.  Let  H[  be  the  part  of  H\  that  is  of  the  first  order  in  the  chosen  expansion 
parameter.  Then,  to  first  order  in  this  parameter,  one  finds 

S'(f,  to)  ~  S' (t,  t0)  -  iS'  (t,  h)  As'  (t',  t0)“1^So(t/,  t0)]  dt' ,  (B5) 

Jto 

from  which,  upon  rotation  back  to  the  original  flavor  basis,  the  evolution  matrix  S(t,to)  is 
obtained.  We  shall  now  consider  the  application  of  the  above  formalism  to  the  cases  when 
the  expansion  parameter  is  either  a  or  s \ 3 . 

B.l  Expansion  up  to  first  order  in  a 

In  this  case,  the  zeroth-order  Hamiltonian  H'0(t)  is  obtained  from  Eq.  (Bl)  by  taking  the 
limit  a  — »  0  (i.e.,  A2i  —> ►  0).  The  zeroth-order  evolution  matrix  S'0(t,to)  can  then  be  written 
as  [26] 

/  u(t,t0)  0  v(t,t0)  e~lScp  \ 

S'(t,t0)=  0  f{t, to)  0  ,  (B6) 

V  -v*(t,t0)elScp  0  u*(t,to)  J 

where  u(t,  to)  and  v(t,  t0)  are  to  be  found  from  the  solution  to  the  two-flavor  problem  governed 
by  the  mass  squared  difference  Am^,  mixing  angle  #13,  and  matter  potential  V(t),  and 

f(t,  to)  =  exp  ^  A  it')  df'j  ,  A  (t)  =  -  +  Vit)  J  .  (B7) 

The  parameters  u  and  v  satisfy  |m|2  +  |u|2  =  1.  The  factor  e“1<Scp  was  pulled  out  from  v  for 
convenience  (this  will  simplify  some  of  the  formulas  given  below). 

The  neutrino  oscillation  probabilities  to  first  order  in  a  in  the  case  of  matter  with  an 
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arbitrary  density  profile  can  be  written  as 


Ppp  =  \u\ 


2  ——  s2 2  Im [uv*  (2  cos  2 0i3 I\  +  sin  20i3  4)] , 

JU 

aA 


(B8) 


aA 


Pe n  =  s23|u|2  +  2  — —  s22  s23  Im[m>*  (2  cos2d13  A  +  sin  20i3 12)\  H — —  sin  20i2  sin  2 023 

.4/ 

x  {cos^cp  Im[n*/*(ci3/3  +  si34)]  +  sin<5Cp  Re[v*  f*(c13I3  +  si3/4)]}  ,  (B9) 


PfiT  =  t  sin2  2023  [1  +  |w|2  -  2R e(uf)] 


aA 


sin2  2023  {  sf2  Im 


+ 


(u  -  f*)v*  ( cos  20i3  A  +  i  sin  20i3  J2 


V 


c12L  —  s12  (  s13 4  +  c13  4  —  sin  20i3  J7 


aA 


-  —  sin  26»i2  sin  26*23 1  cos5Cp  [c23  Im[(w  -  f*)v*(c13I3  +  s13/4)] 

+  Im[(c23|n|2  -  s23  -  cos 2d23  uf)  (c13/9  -  s13  /8)]] 

+  sin  dCp  [c23  Re[(n  -  f*)v*(c13I3  +  s13/4)] 

-  Re[(c23|n|2  -  s23  -  cos2023  uf)  (c13/9  -  s13  J8 
where  the  quantities  J4,  J2, . . . ,  ig  are  given  by 

h  =  f  u*(t',t0)v(t',t0)  dt' ,  4  =  [  [v(t',t0)2  -  u*(t',t0)2}  dt'. 


' to 


I  to 


(BIO) 


4  =  /  u*(t',to)f(t',t0)dt',  h=  v(t',t0)f(t',t0)dt', 


'to 


'to 


h=  \u(t',  t0)  | 2  dt', 


h=  \v(t',t0)\2  dt' , 


'to 


'to 


h  =  \  [ u(t ',  t0)u(t',  to)  +  «*(t',  to)v*(t\  to)]  dt' , 


'*0 


4  =  /  M*(t/,  t0)/*(t',  to)  dt' ,  4  =  /  4t',to)/*(t',to)dt'. 


'io 


't0 


Note  that  the  integrals  I5,  I6,  and  4  are  real,  4  +  4  =  t  — 10,  and  the  integrals  4  and  4  can 
be  obtained  from  I3  and  4,  respectively,  by  substituting  /  — ►  /*. 

In  the  case  of  matter  of  constant  density,  the  parameters  u,  v,  and  /  are  given  by 


Cl3At  .  cos26*13  —  A  .  4i3At 
w(t,  0)  =  cos  — - - h  1 - — - Sill 


C 


13 


.sin 20i3  .  Ci3At 
v(t,  0)  =  — 1 — — - sm  ■ 


C 


(Bll) 

(B12) 
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f(t,  o)  =  exp 


i<1+-4>M 


(B13) 


Inserting  these  expressions  into  the  integrals  ii,  I2, . . . ,  I9  and  the  obtained  results  into  Eqs.  (B8)- 
(B10)  leads  to  Eqs.  (45)-(50). 


B.2  Expansion  up  to  first  order  in  S13 


This  case  was  studied  in  detail  in  Ref.  [28]  (Appendix  A),  which  we  closely  follow  here.  The 
zeroth-order  Hamiltonian  H'0(t)  is  obtained  from  Eq.  (Bl)  by  taking  the  limit  013  — >  0.  The 
zeroth-order  evolution  matrix  S'0(t,t 0)  can  be  written  as 

/  x(t,t0)  y(t,to )  0  \ 

-So  (Mo)  =  I  ~y*(t,t0)  x*(t,t0)  0  ,  (B14) 

\  0  0  g(t,t0)  ) 


where  x(t,  t0 )  and  y(t,  t0)  are  to  be  found  from  the  solution  to  the  two-flavor  problem  governed 
by  the  mass  squared  difference  Am^,  mixing  angle  9 12,  and  matter  potential  V(t),  and 


ff(t,t0)  =  exp 


— i  f  A  (t')  d t!  , 


A(t) 


2  E 


1 

2 


A-mj! 
2  E 


+  V(t) 


(B15) 


The  parameters  x  and  y  satisfy  \x\2  +  | y | 2  =  1. 

The  neutrino  oscillation  probabilities  to  first  order  in  S13  in  the  case  of  matter  with  an 
arbitrary  density  profile  can  be  written  as 

Pee  =  \x  |2,  (B16) 

Pen  =  <4  \y\2  -  sin  2^23  Im(yyC') ,  (B17) 

Pfj.T  =  «23c23  \x  -9* I2  +  sin  2023  Im[(x  -  g*)g(s223B  -  cj3D )] ,  (B18) 


where  the  quantities  B,C,  and  D  are  defined  as  [47] 


B  =  B{t ,  t0)  =  t0)  +  blx>t(t,  t0) ,  (B19) 

C  =  C(t,t0 )  =  a*I**  to{t,  t0)  -  &%%0(t,to) ,  (B20) 

D  =  D(t,  t0)  =  t0)  +  b*I*to(t,  t0) .  (B21) 


Here  the  integral  s  is  defined  as 


0)  =  /  (p(t',s)g{t',s)dt' , 

■'to 


and  the  quantities  a  and  b  are  given  by 

Am2 


31 


2  E 


b  = 


2  E 


(B22) 


Sl3(l-4«)e-Wcp, 

(B23) 

2 

21  — 

— si3si2Ci2e  0CP  . 

(B24) 
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In  the  case  of  matter  of  constant  density,  the  parameters  x,  y,  and  g  are  given  by 


aC  12  At  cos  2  612  — A/ a  aC^At 

x(t,  0)  =  cos - - - 1-  1 - — - Sill  ■ 


c 


12 


.sin 26*12  .  aCV2At 
y{t,  0)  =  — 1 — — - sin 


git,  0)  =  exp 


C\2  L 

,(A  +  a  —  2)  A 


(B25) 

(B26) 

(B27) 


Inserting  these  expressions  into  Eqs.  (B19)  -  (B22)  and  the  obtained  results  into  Eqs.  (B16)- 
(B18)  leads  to  Eqs.  (63)-(68). 
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